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ABSTRACT 

When  a  radar  pulse  impinges  upon  a  target,  the  resultant 
scattering  process  can  be  solved  as  a  linear  time-invariant 
(LTI)  system  problem.  The  system  has  a  transfer  function  with 
poles  and  zeros.  Previous  work  has  shown  that  the  poles  are 
independent  of  the  exciting  waveform  and  target's  aspect,  but 
they  are  dependent  on  the  target's  structure  and  geometry. 
This  thesis  evaluates  the  resonance  estimation  performance  of 
two  signal  processing  techniques:  the  Kumar esan-Tufts 
algorithm  and  the  Cadzow-Solomon  algorithm.  Improvements  are 
made  to  the  Cadzow-Solomon  algorithm.  Both  algorithms  are 
programmed  using  MATLAB.  Test  data  used  to  evaluate  these 
algorithms  includes  synthetic  and  integral  equation  generated 
signals,  with  and  without  additive  noise,  in  addition  to  new 
experimental  scattering  data  from  a  thin  wire,  aluminum 
spheres  and  scale  model  aircraft. 
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I.  INTRODUCTION 

When  a  radar  pulse  impinges  upon  a  target,  the  resultant 
scattering  process  can  be  solved  as  a  linear  time-invariant 
(LTI)  system  problem.  The  system  has  a  transfer  function  with 
poles  and  zeros.  This  description  is  provided  by  the 
Singularity  Expansion  Method  (SEM)  developed  by  Dr.  Carl  Baum 
[Ref.  1],  of  the  Air  Force  Weapons  Laboratory  at  Kirtland  Air 
Force  Base.  Baum's  SEM  describes  the  induced  current  and 
scattered  field  using  damped  sinusoidal  waveforms  which 
resonate  with  complex  natural  frequencies  unique  to  the 
object.  These  frequencies  are  determined  by  the  object's 
composition  and  structural  geometry.  Natural  frequencies  are 
also  the  complex  poles  of  the  transfer  function.  These  poles 
are  independent  of  the  incident  electromagnetic  excitation, 
including  aspect  and  polarization,  as  initially  postulated  by 
Moffatt  and  Mains  [Ref.  2].  Morgan  [Ref.  3]  has  shown  that  a 
target  scattering  response  (following  illumination)  can  be 
represented  as  a  weighted  expansion  of  complex  natural 
resonances  whose  poles  are  independent  of  the  incident 
electromagnetic  excitation.  The  problem  of  classifying  a 
radar  target  can  be  solved,  in  principle,  by  determining  the 
pole  positions  of  the  target's  response. 


Although  this  idea  is  not  new,  early  attempts  to 
demonstrate  the  practicability  of  such  an  identification 
system  have  produced  poor  results  due  to  the  presence  of  noise 
in  the  system.  Two  separate  signal  processing  algorithms  have 
been  developed  by  Kumaresan-Tufts  [Ref.  4]  and  Cadzow- 
Solomon  [Ref.  5]  to  locate  target  poles  with  a  high 
degree  of  accuracy,  in  the  presence  of  noise.  The  Kumaresan- 
Tufts  algorithm  is  intended  for  purely  auto-regressive  (AR) 
signals,  while  the  Cadzow-Solomon  algorithm  is  intended  for 
auto-regressive  moving-average  (ARMA)  signals.  This  thesis 
evaluates  the  viability  of  these  two  signal  processing 
techniques  by  using  new  experimental  electromagnetic 
scattering  data.  The  use  of  this  method  improves 
identification  of  natural  resonances  in  noisy  signals  with  the 
correct  selection  of  problem  parameters. 

A.   THE  PROBLEM 

The  approach  employed  to  classify  radar  targets  is 
comprised  of  two  steps,  based  on  natural  resonances.  The 
first  step  locates  the  position  of  the  poles  of  the  targets  of 
interest.  Numerical  analysis  of  integral  equation  techniques 
will  derive  the  poles  for  simple  targets,  e.g.,  a  thin  wire. 
For  more  complicated  targets  such  as  aircraft,  the  poles  must 
be  extracted  from  actual  measurements  of  the  target's  response 
to  incident  electromagnetic  excitation.    The  information 


collected  can  then  be  used  to  form  a  database  for  comparison 
to  the  data  obtained  in  actual  field  use. 

The  second  step  for  target  classification  is  the 
comparison  of  the  data  obtained  in  field  measurements  with  the 
data  contained  in  the  database.  The  target  is  classified 
based  on  the  closest  data  match.  One  method  of  accomplishing 
the  classification  is  to  use  the  same  signal  processing 
technique  that  was  used  to  form  the  database  initially.  Speed 
is  an  important  factor  in  the  classification  process,  as  the 
time  required  to  achieve  classification  must  be  less  than  the 
time  for  the  target  to  become  a  threat.  Another,  more 
efficient  way  to  perform  the  comparison  is  to  employ  the  use 
of  annihilation  filters,  as  proposed  by  Morgan  and  Dunavin 
[Ref.  6].  This  approach  employs  energy  cancellation 
based  upon  the  location  of  the  target's  poles.  For  example, 
a  total  target  response  (including  early-time  and  late-time) 
may  be  subjected  to  a  bank  of  annihilation  filters.  The 
filter  that  corresponds  to  the  proper  target  will  exhibit,  in 
theory,  a  response  coincident  with  the  driven  portion  of  the 
target  response  while  annihilating  the  signal  during  the  late- 
time  portion  of  the  response. 

A  system  used  for  classification  of  radar  targets  would 
require  a  bank  of  annihilation  filters  corresponding  to  each 
class  of  target.  The  system  could  then  determine  the  class  of 
target,  based  on  the  filter  producing  the  minimal  amount  of 
output  energy. 


B.   BACKGROUND 

Transient  electromagnetic  scattering  can  be  described  by 
showing  an  incident  field  in  free  space  illuminating  a 
perfectly  conductive  finite-sized  object.  Figure  1  depicts 
induced  current  on  a  scattering  body.  The  magnetic  field 
integral  equation  (MFIE) 
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Figure  1     Transient  electromagnetic  scattering 

describes  induced  current  on  the  surface  of  a  scattering  body 
as: 

J(r,t)=2fixHi(r,t)+f[    Id,?,  t)  -J(r,  t-   l7"^  1  )  -35    (1) 

J  J  Spy  C 

where , 

•  J   is  the  electric  current  density, 

•  rf  is  the  unit  normal  vector  outward  to  the  surface, 


•  H1    is  the  incident  magnetic  field  intensity  at  the 
surface, 

•  K   is  the  dyadic  Green's  function  kernel 

and,  the  principal-value  (PV)  type  integral  excludes  the  point 
r^?7 .  The  cross  product,  2i?x#i ,  represents  the  physical 
optics  portion  of  the  induced  current.  The  principal-value 
surface  integral,  Spv,     represents  the  contribution  due  to 

induced  physical  optics  currents  other  than  the  contribution 
from  the  cross  product.  This  current  represents  "feedback" 
effects  from  currents  at  other  points  on  the  object. 

When  ~H1=0 ,  Equation  (1)  solutions  become  the  source-free 

("natural")  current  modes  of  the  scattering  problem.  These 
modes  may  be  described  by  the  sum  of  the  exponentially  damped 
sinusoids,  or:    Jj](r)  exp  (snt)  ,  where  sn    are  the  poles  or 

natural  resonance  frequencies  in  complex  conjugate  pairs  of 
the  form 

Sn=°n+J  '<*n-  <2> 

where , 

•  an   =  the  damping  rate  in  Nepers/sec, 

•  wn  =  the  frequency  in  radians/sec. 

These  resonance  frequencies  are  functions  of  the  geometry  and 
composition  of  the  scattering  object.  Although  the  number  of 
poles  is,  in  theory,  infinite  for  any  given  object,  only  a 


finite  subset  will  be  excited,  due  to  the  finite  bandwidth  of 
the  incident  field. 

Figure  2  represents  a  plane  wave  impulse  illumination, in 
which  part  of  the  object  has  been  illuminated  from  the 
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Figure  2     Plane  wave  impulse  illumination 


incident  field  wavefront.  The  scattered  field  is  composed  of 
two  parts:  the  early-time  driven  response  and  the  late-time 
natural  mode  response.  The  early-time  driven  response  occurs 
as  the  excitation  wavefront  travels  over  the  object  and  each 
point  becomes  excited.  The  late-time  natural  mode  response 
occurs,  by  definition,  when  the  excitation  has  been  completed 
over  the  complete  body.  Throughout  the  early-time  phase,  the 
impulsive  plane  wave  incident  field  will  be  identically  zero 
at  all  points  on  the  surface  except  on  the  conformal  ring 


where  the  intersection  with  the  wavefront  occurs.  This  area 
is  indicated  by  the  dotted  line  in  Figure  2.  This  confonnal 
ring  changes  cross-sectional  shape  and  position  as  the 
wavefront  moves.  The  induced  surface  current  on  the  ring  is 
therefore  composed  of  the  physical  optics  term  plus  the 
feedback  current,  as  described  in  Equation  (1)  .  Thus,  no 
induced  current  is  present  at  the  area  of  the  object  ahead  of 
the  wavefront. 

The  back-scattered  far-field,  resulting  from  the  surface 
current  on  the  object,  will  be  of  the  form 

Hs(-z$,  t)  =— — •-—  [[pxJCr,t  t-   l7"7'!  )  -dS7.  (3) 

4ncr  dtJJ  c 

s 

The  unit  vector  j?  indicates  the  direction  of  the  plane  wave's 
propagation.  The  back-scattered  far-field  for  a  fixed  point 
is  then  a  result  of  integrating  the  current  at  every  point  on 
the  object's  surface.  Thus,  the  back-scattered  far-field  will 
be  of  the  form 


Hs(-r£,  t)=u(t--£)-[Hpo(-r£,  t)  +  £  Hn(-r$,  t)  exp  (snt)  ]  .  (4) 
c  a— « 

n*0 


Two  individual  terms  emerge  from  these  calculations.  The 
first  term  in  Equation  (4)  represents  the  physical  optics 

scattered  field  generated  by  the  2i5x^1  driven  current.   The 

second  term  represents  the  scattered  field  produced  by  the 


source-free  wake  current  or  "feedback"  current.  The  physical 
optics  scattered  field  is  highly  aspect  dependent. 

During  the  early-time  portion  of  the  target* s  response, 
the  scattered  field  produced  by  the  "feedback"  current  term  in 
Equation  (4)  contains  exponential  resonance  terms  with  time- 
varying  coefficients  Hn,     as  described  by  the  Singularity 

Expansion  Method  (SEM)  [Ref .  1] .  This  form  of  the  SEM 
expansion  is  termed  "Class  2"  and  is  presented  analytically  by 
Morgan  [Ref.  7].  For  the  monostatic  scattering  case, 
the  transition  from  early-time  to  late-time  in  the  target's 
response  will  occur  at 

At=T+(l+cosa)-^,  (5) 

c 

after  the  leading  edge  of  the  scattered  pulse  returns  to  the 
radar  antenna,  where, 

•  x  is  the  incident  pulse  width, 

•  D  is  the  target's  largest  dimension,  and 

•  a  is  the  aspect  of  the  target. 

The  term  a  also  represents  the  angle  between  the  target's 
largest  dimension  and  the  direction  of  wave  propagation,  where 
c  is  the  velocity  of  propagation  or  speed  of  light.  At  this 
transition  time,  the  physical  optics  field  vanishes.  During 
the  late-time  portion  of  the  target's  response  only  the  second 
term  in  Equation  (4)  remains.  However,  it  is  now  comprised  of 


8 


constant  coefficients,  HD .   This  form  of  the  late-time  SEM 

expansion  is  termed  "Class  1". 

The  early-time  scattered  field  is  therefore  constituted  of 
both  a  physical  optics  term  and  a  "Class  2"  SEM  expansion  with 
time-varying  coefficients,  while  the  late-time  scattered  field 
is  described  by  a  simple  "Class  1"  SEM  expansion  with  constant 
coefficients. 


II.  POLE  EXTRACTION  ALGORITHMS 

Obtaining  the  target  response  is  the  first  step  to  be 
completed  in  a  Non-Cooperative  Target  Recognition  (NCTR) 
system.  Once  the  response  is  obtained  in  time  domain,  the 
target's  poles  may  be  located.  Numerically,  these  poles  may 
be  described  by  the  damping  rate  a  and  the  radian  frequency 
<i)  or,  in  other  words,  the  complex  number  indicating  the  pole's 
location  in  the  s-plane  or  z-plane.  The  location  of  the  poles 
may  also  be  described  analytically  by  solving  the  boundary- 
value  electromagnetic  problem.  This  may  be  easily  done  only 
for  simple  shape  targets  like  thin  wires  and  spheres. 
However,  for  general  targets,  detailed  knowledge  of  the 
target's  composition  (dimensions  and  materials)  and,  also, 
access  to  a  large  amount  of  computing  power  is  required. 

The  next  section  of  this  thesis  discusses  the  analytical 
method  which  is  used  to  compare  the  poles  extracted  from 
measurements  taken  with  a  thin  wire.  The  viability  of  the 
Cadzow-Solomon  algorithm  will  also  be  tested. 

A.   EARLY  METHODS 

The  NCTR  system  is  used  to  discriminate  and  to  classify 
targets  with  relatively  similar  shapes  and  dimensions.   In 
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such  cases,  the  location  of  the  poles  from  these  targets  will 
be  relatively  close,  requiring  a  high  degree  of  precision  in 
measurement  and  estimation.  A  basic  method  in  signal 
processing  is  to  obtain  the  signal's  spectrum  using  the  Fast 
Fourier  Transform  (FFT) .  The  FFT  algorithm  yields  results  in 
a  short  amount  of  time.  However,  the  FFT  may  be  used  only  as 
a  tool  to  add  information  to  the  basic  methodology  when 
locating  radar  target  natural  resonances.  This  is  because  the 
frequency  resolution,  equal  to  the  reciprocal  of  the  time 
interval  between  the  samples,  is  of  the  order  of  several  MHz 
or  higher.  Another  limitation  of  the  FFT  is  that  it  provides 
only  real  frequencies  within  a  signal,  while  both  the  damping 
rate  and  the  frequency  of  the  poles  are  required  for  target 
classification. 

Thus,  other  methods  providing  more  accurate  results  needed 
to  be  developed. 

1.  Direct  Minimization 

As  previously  described,   the  transient  scattered 
signal  waveform  will  have  the  form 

y ( t)  =yE(  t)  [u( t)  -u(t-T0)  ]  +yL( t)  -u(t-T0)  +N( t)  (6) 

where , 

•  yB(t)    is  the  early-time  portion, 

•  yL(t)    is  the  late-time  portion  and 

•  N(t)    is  the  undesirable  noise  and  clutter. 
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Morgan  [Ref.  7]  presents  the  way  to  determine  the  natural 
resonance  parameters  by  modeling  the  late-time  scattering 
response  as  a  sum  of  damped  sinusoids, 


9L ( t)  =£  AD-exp  (ant)  -cos  («flt+4>n)  (7) 


where  y     represents  the  modeled  signal   and  A^-exp  ( j$n) 

represents  the  aspect  dependent  residues.  The  digital  domain 
form  of  the  model  is 

N 

?(k-A  t)  =?k=Y,  An*exp  (oD-kA  t)  -cos  (o)n-icA  t+4>a)       <«> 

D'l 

where , 

•  AD   is  the  amplitude, 

•  an   is  the  damping  rate, 

•  <*>n   is  the  radial  frequency,  and 

•  <t>n  is  the  phase. 

After  comparing  the  modeled  signal  fn  to  the  actual  received 
discrete  signal  yn,    the  mean  square  error  may  be  obtained  as 

en=(yn~?n)2  •  Tne  four  parameters  of  the  model  (An,  <J>n,  aD,  a>n) 

may  then  be  adjusted  in  order  to  derive  the  mean  square  error 
minimum.  As  this  is  a  multi-dimensional,  highly  non-linear 
minimization  problem,  it  is  computationally  inefficient. 
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2.  Prony's  Method 

Blaricum  and  Mittra  [Ref.  8]  present  a  novel 
approach  for  systematically  deriving  the  complex  poles  and 
residues  of  a  target  from  a  set  of  time  domain  data.  The 
method  is  based  on  Prony's  algorithm,  used  to  model  the  late- 
time  response  of  a  radar  target.  The  set  of  time  domain  data 
is  the  discrete  set  of  sampled  transient  values  of  the  impulse 
response  I(tn)    or  ID.      The  method  is  based  on  the  fact  that 

In   must  satisfy  a  difference  equation  of  order  N,  written  as 

EvJp*— ***-  (9) 

where  p+k=0,l, . . .  ,2N-1  and  2N  are  the  data  samples  used.  This 
leads  to  the  solution  of  a  matrix  equation 


Ix     I2    .. 
I2     J3     .. 


■fw  In+i 


I» 

■ 

N 

*» 

■*W+1 

■ftf+1 

• 

aN-l 

= 

^N+2 

"*2W-1. 

.  *1. 

.     2tf. 

(10) 


After  the  coefficients  ad    are  found,  the  method  solves  the 
roots  of  a  polynomial  as 


z^z*-1- 


a„=0. 


(11) 


If  dm  (m=l,2/ . . . ,N)  represents  those  roots,  then  the  poles  of 
the  system  model  in  the  z-plane  are  the  roots  dm.  The  poles 
in  the  s-plane  may  be  found  using  the  formula 
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v-^-.  <»> 


where  A  t  is  the  time-stepping  interval  used  in  obtaining  the 
sampled  data.  The  matrix  in  Equation  (10)  is  in  the  form  of 
a  transposed  Vandermonde  matrix,  whose  inverse  can  be  computed 
in  closed  form.  This  method  requires  at  least  2N  equally 
spaced  transient  data  samples  to  find  N  poles.  If  greater 
than  2N  samples  are  desired,  then  one  may  obtain  a  least 
squares  type  fit  to  the  matrix  in  Equation  (10) .  Since  the 
system  order  is  not  known  a  priori  in  the  NCTR  problem, 
Blaricum  and  Mittra  [Ref.  9]  present  two  schemes  for 
determining  the  number  of  poles  contained  in  the  transient 
response.  The  first  is  the  Householder  orthogonal ization 
method,  and  the  second  is  an  eigenvalue  method.  When  the  data 
are  noisy,  Blaricum  and  Mittra  [Ref.  9]  indicate  that  the 
eigenvalue  method  is  the  better  method.  This  method  will  be 
discussed  in  Section  3  of  this  chapter  when  a  bias 
compensation  method  is  examined. 
3.  Kumaresan-Tufts  Algorithm 

Kumaresan  and  Tufts  [Ref.  4]  modified  Prony's  method 
by  using  the  "backward  linear  prediction"  technique  and 
"singular  value  decomposition"  (SVD)  to  alleviate  the 
sensitivity  to  noise  and  the  need  for  a  priori  knowledge  of 
the  system  order. 
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Analytically,  Kumaresan-Tufts  algorithm  modifies 
Prony's  method  as  follows: 

•  The  system  order  is  intentionally  overestimated.  The 
excess  order  provides  the  flexibility  to  the  system  to 
model  the  noise,  improving  the  estimation  of  parameters  of 
exponentially  damped  sinusoidal  signals  in  noise. 

•  Singular  value  decomposition  alleviates  severe  ill- 
conditioning  of  the  data  matrix. 

•  The  causality  of  the  system  is  used  to  separate  the 
computed  signal's  poles  from  the  spurious  computed  noise 
poles  introduced  by  the  overestimated  system  order. 

The  separation  of  the  signal  poles  from  the  noise 
poles  is  the  result  of  the  "backward  prediction"  that  causes 
the  signal  poles  to  fall  outside  the  unit  circle  in  the  z- 
plane,  while  the  extraneous  "noise  poles"  fall  inside  the  unit 
circle. 

Moderately  large  values  of  system  order  are  essential 
in  improving  the  accuracy  of  the  pole  location  estimates. 

a.  System  Model 

Suppose  that  the  N  samples  of  the  observed  time 
domain  data  of  a  response  signal  y(n)  consists  of  samples  of 
M  exponentially  damped  signals  in  white  Gaussian  noise  w(n)  , 
as  represented  by 

M 

y(n)  =J2  a^ex.piSffl)  +w{n)  ,  n=0,l,  .  .  .  ,N-l.  (13) 

The  following  linear  prediction  equations  can  be  formed,  using 
real  time  domain  data  in  the  backward  direction. 
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y(2) 
y(3) 

y(N-L+l) 


y(3) 
y(4) 


yU+D 
yU+2) 

y(tf) 


a(D 
a  (2) 

a(L) 


y(D 
y(2) 

y(isr-L) 


(14) 


The  terms  of  the  above  equation  may  be  also  represented  by 

D-a=-h  (14a) 

where , 

•  D   is  the  data  matrix  (N-L) xL, 

•  a  is  the  vector  of  backward  prediction  coefficients 
(Lxi) ,  and 

•  h   the  data  vector  (N-L)xi. 

In  the  above  matrix  equation,  L  represents  the  overestimated 
system  order,  chosen  to  satisfy  the  inequality  M<L<N-M.  The 
matrix  equation  is  solved  with  the  SVD  method,  using  the 

pseudoinverse  matrix  D* ,  as  a=D*-(-h)   .  Here  the  coefficients 

a±   correspond  to  those  in  Equation  (11) ,  where  the  roots  of 

the  polynomial  define  the  poles  of  the  system  model  in  the  z- 
plane. 

As  with  Prony's  method,  the  Kumaresan-Tufts  method 
is  intended  for  purely  auto-regressive  (AR)  signals.  Both 
methods,  therefore,  use  the  late-time  portion  of  the  response 
signal  when  addressing  the  target  classification  problem. 
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b.  Singular  Value  Decomposition 

Singular  value  decomposition  (SVD)  is  the  basic 
technique  on  which  the  Kumaresan-Tufts  algorithm  is  based. 
There  are  two  main  advantages  with  SVD: 

•  It  helps  to  alleviate  the  effects  of  ill-conditioning  of 
the  data  matrix. 

•  It  separates  the  signal  poles  from  the  noise  poles. 

The  following  discussion  of  SVD  is  a  synopsis  of  material 
taken  from  Strang  [Ref.  10]. 

The  SVD  is  closely  associated  with  the  eigenvalue- 
eigenvector  factorization  of  a  symmetric  matrix  :  A=Q-A-QT. 
SVD  makes  a  similar  factorization,  but  for  any  (mxn)  matrix 
A,    as  follows: 

a=qx*-q2t  (is) 

where  Qx,    Q2   are  two  orthogonal  matrices  and  S  is  the  diagonal 

(but  rectangular)  matrix  with  its  positive  entries  (also 
called  sigma,  in  the  form  of  ax,  o2,  .  .  . ,  az) .  These  entries  are 

the  singular  values  of  A,  filling  the  first  r  places  on  the 
main  diagonal  of  E,  where  r  is  the  rank  of  A.  These  entries 
sigma  are  also  the  square  roots  of  the  nonzero  eigenvalues  of 
both  AA  T   and  ATA .  The  columns  of  Qx    (mxm)  are  eigenvectors  of 

AAT ,    and  the  columns  of  Q2    (nxn)  are  eigenvectors  of  ATA. 
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The   SVD   works   well   for   numerically   stable 
computations.    An  initial  reason  is  that  Qx     and  Q2     are 

orthogonal  matrices  that  never  change  the  length  of  a  vector. 
Since  ||Cx||2=xr£>r£x=||.xj|2 ,  multiplication  by  Q  cannot  destroy 
the  scaling.  A  second  reason  is  that  SVD  gives  a  more  stable 
measure  of  the  rank  of  A. 

The  prime  use  of  SVD  is  to  solve  every  linear 
system  in  the  form  of  Ax=b .  For  every  matrix  A    (m  x  n)  ,  which 

can  be  factored  into  A=Q1*L'Q2T ,  another  matrix  can  be  defined. 
This  matrix  will  be  the  pseudoinverse  of  A,  as  follows: 

A+=Q2WiT  <16> 

where,  Qx,  Q2   are  the  same  orthogonal  matrices  found  with  the 

SVD.   The  value  S*  is  the  diagonal  (n  x  m)  matrix,  with  its 

11       1 
entries  on  the  main  diagonal  being   — ,  — ,..., — .    The 

°l  °2  °r 

optimal  solution  of  Ax=b   is 

x*=A*-b  (17) 

that  is,  the  minimum  length  solution  to  the  nearest  equation, 
Ax=p.  As  the  column  space  of  A*  is  in  the  row  space  of  A, 
then  x*   is  always  in  the  row  space  of  A . 
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To  solve  the  system  of  equations  in  Equation  (14), 
the  pseudoinverse  of  matrix  D  may  be  found  as  in  Equation  (16) 
in  the  form  D+  =  W*'UT ,  where  2+  will  be  a  (L  x  (N-L) )  matrix, 
whose  L  singular  values  are  the  reciprocals  of  those  found  in 
the  2  matrix.  The  minimum  length  least  squares  solution  of 
D-a=-h  is  a=D+-{-h)=VWT-(-h)  . 
c.  Bias  Compensation 

In  the  discussion  of  the  SVD  above,  S  is  a  (N-L)  xL 
diagonal  matrix,  with  its  entries  being  the  square  roots  of 
the  nonzero  eigenvalues  of  both  DDT  and  DTD.  In  the  noiseless 
case,  the  prediction  equations  are  satisfied  exactly.  If  the 
system  order  has  been  overestimated  as  L,  when  M  is  the  actual 
order  of  the  system,  the  diagonal  of  S  splits  into  a  signal 
subspace  with  M  positive  singular  values  of  <s1,o2, . . . ,oM   and 

a  subspace  with  L-M  zero  values.  For  the  noisy  signal  case, 
the  first  M  positive  signal  values  are  perturbed  into  a  noisy 
signal  subspace 

olfa2,  ...  ,aM  with  eigenvalues  X^X^.  .  .  z.X'M ,  (18) 

and  a  noise  subspace  with  L-M  values 

°M*i'au*2'  •  '  •  •  °l  with  eigenvalues  \'M+1z\'M+2z •  •  •  *\'L  .      (19) 

Kumaresan  and  Tufts  noticed  [Ref .  4]  that  the  noise 
perturbed  the  signal's  singular  values,  the  reason  for  the 
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bias  towards  the  unit  circle  in  the  z-plane  in  the  pole 
position  estimates.  Kumaresan  and  Tufts  proposed  a 
compensation  method  which  moves  the  poles  back  towards  the 
center  of  the  z-plane.  This  method  is  based  on  averaging  the 
smallest  L-M  singular  values,  °M*i'aM+2'  •  •  • » °l»  and  then 
subtracting  this  average  value  from  each  of  the  first  M 
singular  values,  alta2, . . . ,aM.       The  smallest  L-M  values  are 

then  set  to  zero  and  a  new  matrix,  2  ,  is  used  to  compute  the 
pseudoinverse  D* .  Although  Kumaresan  and  Tufts  claim  that 
this  method  gives  better  results,  the  testing  completed  for 
this  thesis  shows  that  the  bias  towards  the  center  of  the  z- 
plane  is  greater  than  the  bias  error  that  the  noise  makes. 
This  result  causes  much  more  perturbation  on  the  pole  position 
estimates. 

Norton  [Ref.  11]  proposed  another 
compensation  method  based  on  the  Eigenvalue  Shifting  Theorem. 
The  noisy  data  matrix  may  be  described  by  D=S+N,  where  S  is 
the  noiseless  signal  data  matrix  and  N  is  the  noise  data 
matrix  of  the  wide-sense  stationary  white  noise  process,  v^^ , 

written  as 


N- 


VN-L    -    VN. 


(20) 
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The  expected  value  of  DDT   may  be  found  from 

DDT=E[(S+N)  (S+N)*]  =E[SST]  +E[SNT]  +E[tfSr]  +E[NNT]  .      (21) 

As  white  noise  has  a  mean  of  zero,  the  two  terms  EiSN*]  and 
EiNSF]  become  zero.  Since  the  noise  is  wide-sense  stationary, 
ElNN7]  =al'I ,  where  ov  is  the  noise  variance  and  J  is  the 

identity  matrix.  As  S  is  deterministic,  EiSS7]  =SST.  This 
leads  to  the  formula 

E[DDT]=SST+ol-I.  (22) 

The  eigenvalue  shifting  theorem  [Ref.  11]  describes  the  case 
when  the  eigenvalues  of  SS T   are  X± ,   where  the  eigenvalues  of 

EiDD7*]    are  A^+oJ .   Therefore,  the  eigenvalues  of  DDT,    which 

are  the  squares  of  the  singular  values  of  D,  are  increased  by 
the  noise  variance,  oj . 

In  the  noiseless  case,  these  singular  values  would 
be  zero.    If  the  noise  singular  values  are  squared  and 

averaged,  an  unbiased  estimation  of  Oy  may  be  obtained. 

Norton  [Ref.  11]  proposed  correcting  the  signal  singular 
values  by  first  subtracting  the  noise  variance  al    from  the 

squares  of  the  first  M  diagonal  entries  of  matrix  S  and  then 
taking  the  square  root  of  the  reduced  values  as  the  new 
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diagonal  entries  of  matrix  S  (while  the  noise  singular  values 
are  set  to  zero)  .  As  in  the  Kumaresan-Tufts  bias  compensation 
method,  the  pseudoinverse  D*   can  be  found  from  the  new  matrix 

S. 

Although  Norton's  method  seems  to  be  a  more  correct 
bias  compensation  method,  both  methods  are  based  on  the  fact 
that  a  decision  has  to  be  made  about  the  actual  order  of  the 
system.  The  separation  between  the  signal  eigenvalues  and  the 
noise  eigenvalues  is  not  readily  obvious,  because  the 
eigenvalues  of  matrix  DDT  or  DTD  appear  to  steadily  decrease. 
d.  Performance  and  Earlier  Results 

Kumaresan  and  Tufts  tested  the  algorithm,  obtaining 
good  results  using  synthetic  data  with  various  levels  of  white 
Gaussian  noise,  as  low  as  15  dB  SNR.  Norton  [Ref.  11] 
developed  the  algorithm  as  a  computer  subroutine  and  tested  it 
with  various  sets  and  types  of  data.  Synthetic  data  were 
tested  at  various  SNR  values,  ranging  from  90  dB  to  7  dB. 
When  the  SNR  decreased,  the  error  in  pole  position  estimates 
increased.  Norton  claimed  that  the  algorithm  gave  good 
results  for  SNR>20  dB.  Norton  also  used  simulations  of  thin 
wires  produced  by  Morgan's  Time-Domain  thin  wire  integral 
equation  (TDIE)  computer  program  to  test  the  algorithm.  The 
results  of  those  tests  illustrated  the  aspect  independence  on 
estimated  poles  of  the  target,  [Ref.  11] . 
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Larison  [Ref.  12]  programmed  the  algorithm 
using  Fortran  and  tested  synthetic  and  thin  wire  integral 
equation  data  at  various  SNR's  ranging  from  90  dB  to  7  dB. 
Like  Norton,  he  claimed  good  results  in  extracting  the  low 
frequency  poles  position  with  SNR  as  low  as  2  0  dB  by  using 
Norton's  bias  compensation  method.  Larison  maintained  that 
the  two  most  critical  parameters  are: 

•  to  select  the  appropriate  starting  point  to  begin  the  data 
processing 

•  to  select  the  appropriate  system  order  to  provide  the  best 
possible  results. 

The  algorithm  has  been  tested  by  the  author  of  this 
thesis,  after  developing  the  algorithm  using  MATLAB.  The 
algorithm  was  tested  using  synthetic  data  at  various  SNR's 
ranging  from  30  dB  to  10  dB,  without  using  the  bias 
compensation  method.  The  synthetic  signal  response  is  based 
on  ten  pole  pairs  within  a  frequency  range  of  1-10  GHz,  with 
a  medium  Q  damping  factor  using  k=0.7.  This  chosen  value  of 
k  simulates  typical  expected  levels  of  damping  from  measured 
scattering  responses  of  real  targets  (e.g.,  the  thin  wire). 
Damping  factor  and  SNR  level  are  discussed,  in  detail,  in 
Chapter  III,  Section  A,  of  this  thesis.  Table  I  shows  the 
poles  and  residues  used  for  the  testing. 

Figures  3  to  10  illustrate  the  evaluation  of  the 
Kumaresan-Tufts  algorithm.  Figures  3  to  6  illustrate  the 
synthetic  signal  waveform  for  various  SNR's.  Figures  7  and  8 
illustrate  the  spectrum  of  the  synthetic  signal  for  SNR=30  dB 
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and  SNR=10  dB.  Figures  9  and  10  illustrate  the  true  and  the 
algorithm  poles  in  both  the  s-plane  and  z-plane,  for  various 
SNR's.  For  the  above  simulation,  the  values  used  for  N  and  L 
were  300  and  60,  respectively,  while  M  (actual  order)  was  20. 
The  simulation  revealed  that  as  the  noise  increases  relative 
to  the  signal,  there  is  a  bias  of  all  the  poles  towards  the 
unit  circle,  specially  the  higher  frequency  poles. 

Table  I  MEDIUM  Q  SYNTHETIC  POLES 


fn 
[GHz] 

[GNp/s] 

0) 
n 

[Grad/s] 

An 

0 

1 

-0.3562 

6.2752 

1 

0 

2 

-0.7124 

12.5504 

1 

0 

3 

-1.0687 

18.8256 

1 

0 

4 

-1.4249 

25.1007 

1 

0 

5 

-1.7811 

31.3759 

1 

0 

6 

-2.1373 

37.6511 

1 

0 

7 

-2.4935 

43.9263 

1 

0 

8 

-2.8498 

50.2015 

1 

0 

9 

-3.2060 

56.4767 

1 

0 

10 

-3.5622 

62.7518 

1 

0 
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Figure   3  Generated  synthetic  signal  without  noise 
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Figure   4  Generated  synthetic  signal  with   30  dB  SNR 
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Figure   5  Generated  synthetic  signal  with  20  dB  SNR 
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Figure  6     Generated  synthetic  signal  with  10  dB  SNR 
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Spectrum  of  the  SNR=30  dB  synthetic  signal 
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Figure  8  Spectrum  of  the  SNR=10  dB  synthetic  signal 
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Figure   9  Kumaresan-Tufts  poles   in  the   s-plane 
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Figure   10         Kumaresan-Tufts  poles   in  the   z-plane 
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B.   CADZOW- SOLOMON  ALGORITHM 

Many  investigators  have  used  both  excitation  and  response 
data  to  identify  linear  systems.  Cadzow  and  Solomon  [Ref .  5] 
proposed  an  identification  procedure  in  which  both  the 
excitation  and  response  are  contaminated  by  white  noise.  This 
method  is  based  upon  the  null  space  characterization  of  an 
associated  "data  matrix" .  The  Cadzow-Solomon  algorithm  is 
used  in  this  thesis  to  simulate  and  demonstrate  the 
identification  problem. 

1.  System  Model 

The  Cadzow-Solomon  algorithm  is  based  upon  the  Auto- 
Regressive  Moving-Average  (ARMA)  model.  In  an  ideal  modeling 
situation,  the  assumption  may  be  made  that  the  excitation 
signal  x(n)  and  the  response  signal  y(n)  are  perfectly 
related  by  means  of  an  ARMA  relationship,  such  that 


L  K 

^aiy(n-i)+y 


y(n)  =V  aiy(n-i)  +£  bjx(n-i) ,  (23) 


as  associated  with  the  transfer  function 


H(z)=     °   *  ■ ^—  (24) 


l-a^-1-.  .  .  -aLz~L 


where , 


a±   are  the  coefficients  which  correspond  to  the  poles  of 
the  transfer  function,  and 
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•  jbj  are  the  coefficients  which  correspond  to  the  zeros  of 
the  transfer  function. 

The  order  of  the  numerator  of  the  system's  transfer  function 

is  K   and  the  order  of  the  denominator  is  L.       The  classic 

modeling  problem  is  to  identify  the  system's  coefficients,  a± 

and  bif    from  a  finite  set  of  observations  of  the  excitation, 

x(n)  ,  and  response,  y(n)  ,  time  series. 

Cadzow-Solomon  present  the  set  of  equations  for  the 
algorithm  in  matrix  form  as  follows: 


y(0)  y(-l)  ...  y(-p) 
y(l)      y(0)   ...  y(l-p) 

y(N)   y(N-l)   ...  y(N-p) 


x(0)    x(-l)  ...  x(-q) 
x(l)      x(0)      ...  x(l-q) 

x(N)    x(N-l)    ...  x{N-q). 


(25) 


or,  more  concisely 


W*A 


(26) 


where , 

•  Xg   is  the  (N+l)x(q+i)  "excitation  data  matrix", 

•  Yp   is  the  (N+l)x(p+i)  "response  data  matrix", 

•  ap   is  the  (p+l)xi  "auto-regressive  parameter  vector"  and 

•  bg   is  the  (q+l)xi  "moving-average  parameter  vector". 
The   least   squares   solution   for  a   causal,   real   data, 
overestimated  ARMA  model  (pzip,  q^q)    can  be  found  by  taking  a 
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linear  algebraic  approach  based  on  the  eigenvalue-eigenvector 
decomposition  of  the  data  matrix  Dp  g=  [Yp:.  -Xg]    as 


DP^'Dp,q'uk=Xk'uk,   lsicsp+g+2. 


(27) 


The  parameter  vector  is  found  from  the  procedure  to  be 


L  «J 


pg2  kd)  I2 


Jr-l 


-1 


****     Uk(l)    _ 


** 


(28) 


where  uk   is  the  eigenvector  associated  with  the  eigenvalue  Xk. 

If  the  actual  system  order  is  (p,  g) ,  then  the  system  order  is 
less   than   the   overestimated   order   (p,  g)  .  At   least 

s=l +min (p-pt  g-g)  of  the  eigenvalues  in  the  decomposition,  as 
described  in  Equation  (27)  ,  must  be  identically  zero  for  the 
noise-free  case. 

Following  the  technique  of  "backward  prediction"  as 
used  in  the  Kumaresan-Tufts  algorithm,  the  matrix  equation 
(Equation  25)  may  be  modified  as  follows: 


y(2) 

..  y(L+D 

x(i)     . 

..     x{K+l) 

y(3) 

..  y(L+2) 

x{2)      . 

..      x(K+2) 

y(N-L+l)    . 

..   ym 

x(N-L)    . 

..  x(N-L+K) 

yd) 
y(2) 

y(N-L) 


(29) 


or,  in  matrix  notation 


uw 


=y. 


(30) 
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This  matrix  is  associated  with  the  coefficients  a±   and  b±   as 

shown  in  Equations  (23)  and  (24) .  The  size  of  the  excitation 
and  response  data  matrices  are  (N-L)x(K+l)  and  (N-L) xL, 
respectively,  while  the  size  of  the  ARMA  parameter  vector  is 
(L+K+l)xl. 

As  in  the  Kumaresan-Tufts  case,  singular  value 
decomposition  can  also  be  used  in  Cadzow-Solomon's  algorithm 
to  provide  the  minimum-norm  solution.  Its  use  reduces  ill- 
conditioning  of  the  data  matrix,  D^,   and  separates  the  signal 

poles  from  the  noise  poles  across  the  unit  circle.  The 
optimal  solution  of  Equation  (30)  is 


■PWV.  (31) 


2.  Bias  Compensation 

When  the  system's  order,  (p,  g)  ,  exceeds  the  true  order 
of  a  noise  contaminated  system,  (p,  g)  ,  for  p>p  and  g>g, 
Cadzow  and  Solomon  maintain  that  there  will  be 

s=l+min(p-p,g-g)  f32) 

eigenvalues,  which  will  asymptotically  approach  (N+l)  •  a2w   for 
large  N,  where  o2v   is  the  noise  variance.   The  eigenvalues  in 
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the  noise-free  case  are  zero,  and  therefore,  the  parameter 
vector  may  be  found  from  the  equation 


(^)bt 


ElMu 


*«o 


•EU*<1)T^ 


(33) 


Jc-0 


where  the  uk   terms  are  the  eigenvectors  associated  with  the 
smallest  (s)  eigenvalues,  Xk. 

However,  the  above  method  may  not  be  applied  exactly 
as  described  above  because 

•  the  excitation  noise,  av,  is  zero  for  the  radar  target 
classification  problem,  and 

•  the  Q-Q,  cannot  be  directly  determined  since  the 
extraneous  zeros  cannot  be  separated  from  the  signal  zeros 
in  the  same  manner  as  the  poles  are. 

Another  problem  occuring  in  the  Cadzow-Solomon's 

algorithm  is  that  additive  noise  is  different  for  input  and 

output  data.   Norton  [Ref.  11]  described  that  if  the  input 

data  noise  is  v±   and  the  output  data  noise  is  wi ,  the  noisy 

data  matrix  may  be  expressed  as 


DyxiPyiDz]*Syx+Nyx 


(34) 


where ,    Nyx=[Ny  \  JVX] ,    or  in  matrix  form 


**= 


vi     •••     Vjfl 


VN-L    -    VN+K-L 


(35) 
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Ny= 


W2         ■"  WL+1 


WH-L*1    "'       WN 


(36) 

Therefore,  the  correlation  product  will  be 

^D^S^S^N^N^  (37) 

As  the  additive  noise  is  different  for  the  input  and  output 
data,  the  noise  correlation  matrix  ElN^N^]    is  not  of  the  form 

a2 1.  So,  the  eigenvalue  compensation  method  that  Norton 
described  as  being  applicable  in  the  Kumaresan-Tufts  algorithm 
is  not  applicable  in  the  Cadzow-Solomon  algorithm. 

By  examining  the  diagonal  entries  of  the  data  matrix 

product,  Dpig'Dp  g,    it  may  be  seen  that  the  first  (g)  entries 

are  close  to  1.  In  the  noise-free  case,  these  first  (g) 
diagonal  entries  are  exactly  equal  to  1.   This  is  a  direct 

consequence  of  the  fact  that  the  diagonal  entries  of  DpTig'Dp  q 
contain  bias  errors  proportional  to  the  noise  variances  o2v  and 

0*.    Since  the  excitation  data  used  in  radar  target 

classification  problems  are  noise-free,  these  bias  errors  are 
proportional  only  to  o2w.         The  bias  compensation  may  be 

obtained  by  setting  the  first  (g)  diagonal  entries  equal  to  1 
and  then  finding  the  eigenvalues  and  eigenvectors  from  the 
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corrected  matrix,  Dp>g'Dp  g.      The  parameter  vector  may  then  be 
found  from 


aP 


.*« 


Pg2  1^(D|2 

*«1      ** 


"Wa  Ujt(l)_ 


E  -nr1"*  <38> 


3.  Earlier  Results 

Norton  [Ref.  11]  programmed  the  Cadzow-Solomon 
algorithm  in  Pascal  and  tested  it  on  various  types  of  data. 
Using  synthetic  data,  Norton  tested  the  algorithm  using  three 
SNR  values,  30  dB,  20  dB,  and  10  dB.  The  results  were  good  at 
the  30  dB  level,  but  as  the  SNR  decreased,  the  error  in  the 
pole  position  estimates  increased.  Norton  also  used 
simulations  of  thin  wire  scattering  produced  by  Morgan's  TDIE 
computer  program  to  test  the  algorithm,  obtaining  better 
results  than  those  obtained  from  the  Kumaresan-Tufts 
algorithm.  These  results  occured  because  the  Cadzow-Solomon 
algorithm  takes  into  account  only  the  input  and  output  data 
and  does  not  require  purely  late-time  data. 

Larison  [Ref.  12]  programmed  the  Cadzow-Solomon 
algorithm  in  Fortran  and  tested  it  using  synthetic  data  and 
TDIE  data,  at  various  SNR's,  ranging  from  90  dB  to  7  dB. 
Using  bias  compensation,  Larison  claimed  good  results  in 
extracting  the  low  frequency  pole  positions  with  SNR  as  low  as 
20  dB. 
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Murphy  [Ref.  13]  tested  the  Cadzow-Solomon 
algorithm  using  Larison's  Fortran  programs,  using  synthetic 
data  as  well  as  thin  wire  integral  equation  generated  data  and 
measured  data.  Murphy  generated  three  separate  signals,  each 
based  on  ten  pole  pairs  covering  the  frequency  range  of  1-10 
GHz.  Depending  on  the  level  of  damping  for  each  signal, 
Murphy  obtained  an  average  error  between  the  true  and  the 
extracted  poles  having  values  on  the  order  of  10"2  for  the 
values 

•  SNR=10  dB  and  for  HIGH  Q, 

•  SNR=2  0  dB  and  for  MEDIUM  Q,  and 

•  SNR=3  0  dB  and  for  LOW  Q. 
Murphy  made  the  observations  that, 

•  The  best  results  were  obtained  by  choosing  a  starting 
point  located  within  several  points  of  the  zero  crossing 
nearest  to  the  first  obvious  response  of  the  excitation. 

•  The  best  results  were  obtained  by  using  a  data  matrix  in 
which  the  overestimated  number  of  poles  to  the  true  number 
of  poles  was  of  the  order  of  2.5,  while  the  number  of 
asking  zeros  was  equal  to  the  number  of  true  zeros. 

•  The  best  results  were  not  always  obtained  when  using  the 
bias  compensation  scheme  as  proposed  by  Norton. 

When  processing  the  thin  wire  data,  Murphy  calculated  the 

feed-forward  order  of  the  system  by  determining  the  length  of 

early-time  as  equal  to  2L/c.  Recent  research  at  NPS,  that  is 

as  yet  unpublished,  indicates  that  this  method  is  not  correct, 

as  each  excited  point  along  the  target  excites  its  adjacent 

point.   This  early-time  may  be  described  as  (1+cosa) * L/c,  as 

shown  in  Chapter  III,  where  a  is  the  aspect  of  the  target. 
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When  processing  the  TDIE  data,  Murphy  did  not  obtain  good 
results  for  SNR=20  dB,  in  spite  of  the  fact  that  the  poles 
extracted  from  thin  wire  measured  responses  appeared  to  give 
good  results  for  the  low  frequency  poles. 
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III.  ALGORITHM  TESTING 

The  initial  objective  of  this  thesis  was  to  evaluate  the 
viability  of  the  Kumaresan-Tufts  and  the  Cadzow-Solomon 
algorithms.  In  a  radar  target  classification  problem  the 
finite  duration  excitation  signal  produces  both  early-time  and 
late-time  response  signals.  In  this  situation,  the  Cadzow- 
Solomon  algorithm  is  more  significant.  Thus,  the  main  effort 
of  this  thesis  was  changed  to  evaluate  and  improve  the  Cadzow- 
Solomon  algorithm  using  both  synthetic  and  real  data.  The 
Kumaresan-Tufts  algorithm  was  evaluated  and  tested  only  with 
synthetic  data,  as  presented  in  Chapter  II. 

The  Cadzow-Solomon  algorithm  was  tested  in  two  phases. 
The  first  phase  of  testing  used  synthetic  data,  while  the 
second  phase  was  performed  with  thin  wire  measurement  data. 
The  synthetic  data  testing  phase  attempted  to  simulate  the 
conditions  expected  from  the  response  of  a  simple  target 
during  the  presence  of  a  stationary  white  noise.  The  thin 
wire  data  testing  phase  attempted  to  evaluate  the  conditions 
appearing  from  a  real  target  response.  Those  conditions  were 
then  compared  with  those  simulated  from  the  computed  Time 
Domain  Integral  Equation  (TDIE)  thin  wire  response. 
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A.   SYNTHETIC  SIGNAL  MODEL 

As  the  Cadzow-Solomon  algorithm  is  based  on  the  ARMA 
model,  the  representation  in  Equation  (23)  has  been  used  to 
produce  the  synthetic  signal  response. 

This  signal  response  is  based  on  ten  pole  pairs  within  a 
frequency  range  of  1-10  GHz,  with  a  medium  Q  damping  factor 
using  k=0.7.  The  poles  were  developed  in  accordance  with  the 
following  equation: 


G)„   271 


(39) 


Table  II  lists  the  s-plane  poles  used  in  the  synthetic  signal 
Table  II   MEDIUM  Q  SYNTHETIC  POLES 


[GHz] 

[GNp/s] 

[Grad/s] 

1 

-0.3562 

6.2752 

2 

-0.7124 

12.5504 

3 

-1.0687 

18.8256 

4 

-1.4249 

25.1007 

5 

-1.7811 

31.3759 

6 

-2.1373 

37.6511 

7 

-2.4935 

43.9263 

8 

-2.8498 

50.2015 

9 

-3.2060 

56.4767 

10 

-3.5622 

62.7518 

The  chosen  value  of  k  simulates  typical  expected  levels  of 
damping  from  measured  scattering  responses  of  real  targets, 
(e.g.,  the  thin  wire  and  scale  model  aircraft  targets).  The 
sampling  frequency  used  to  convert  the  s-plane  poles  to  z- 
plane  poles  was  51  GHz,  based  on  N=256  samples  over  a  time 
window  of  t0=5  nsec: 
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f  _    1  .  N-l  .... 

f'-It~T  (40> 


1.  Coefficient  Generator  and  Recursive  Signal  Generator 

The  coefficients  a±    are  obtained  by  multiplying  the 

terms  {z-zj  (z~z2)  .  .  .  (z-z2Q)  ,  where  zi  is  the  ith  z-plane  pole 
associated  with  the  s-plane  pole  in  the  relationship 

zi=exp[  (Oj+jc^)  -At]  ,  (41) 

and  equal  to  the  product  of  those  terms  with  the  polynomial 

z^-a^-a^16-.  .  .-a20.  (42) 

The  coefficients  b±     are  obtained  by  the  inverse  partial 

R 

fraction  expansion  using  the  term  — ,    where  R±    is  the  ith 

z-z± 

residue  of  amplitude  value  1  and  phase  difference  0  for  each 
of  the  signal  poles.  The  time  domain  signal  response  of  the 
ARMA  model  is  generated  via  Equation  (23)  with  the  values  L=20 
and  K=19.  This  procedure  has  been  developed  in  the  MATLAB 
computer  program  ARMAGEN1.M.  The  program  listing  appears  in 
Appendix  A. 
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2.  Double  Gaussian  Smoothing  Function  Generator 

The  excitation  signal  chosen  for  generating  the  test 
signal  response  was  the  double  Gaussian  waveform,  via  the 
equation 


x(n)  =A1exp(-a1t2)  -A2exp(-a2t2)  , 


(43) 


with 


41n(10) 


Tl 


7^=0.147  nsec 


(44) 


<32  = 


_ 4ln(10) 


T2=0  .314  nsec 


(45) 


where , 


^i= 


aX~^2 


>        **2  ~     1 


(46) 


This  waveform  is  a  wide  Gaussian  pulse  with  a  ten  percent 
width  of  T2   nsec  subtracted  from  a  narrow  Gaussian  pulse  with 

a  ten  percent  width  of  Tx   nsec.   This  results  in  a  bandwidth 

of  1  to  10  GHz.  Figure  12  illustrates  the.  double  Gaussian 
waveform  and  Figure  13  illustrates  the  spectrum  of  the 
waveform.  This  procedure  has  been  developed  in  the  computer 
program  EXCGEN.M  in  MATLAB.  The  program  listing  appears  in 
Appendix  B. 
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3.  Synthetic  Noise  Generator 

Synthetic  noise  was  generated  to  contaminate  the 
signal  response  by  adding  a  time  series  noise  signal  to  the 
signal  response.  The  noise  was  assumed  to  be  wide  sense 
stationary  and  white.  To  produce  a  Gaussian  distribution,  a 
normal  distribution  function  was  multiplied  with  a  standard 
deviation  value  o ,  computed  via  the  equation 

£  [y(k)]2 

o2=  variance  =  i=i • ^-  (47) 

N  SNR 

10  10 

where, 

•  y(n)    is  the  signal  response  data, 

•  N   is  the  number  of  time  points  and 

•  SNR    is  the  Signal-to-Noise  Ratio  in  dB. 

This  procedure  was  developed  in  the  computer  program 
NOISEGEN.M  in  MATLAB.  The  program  listing  appears  in 
Appendix  C. 

4.  Spectrum  Estimation 

Estimation  of  the  power  spectral  density  (PSD) ,  also 
called  the  spectrum  of  the  sampled  signal  response,  is 
obtained  employing  the  Fast  Fourier  Transform  (FFT) .  This 
spectrum  may  be  computed  via  the  equation 

S(k)=±-\Y(k)\2  (48) 

N 
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where , 

•  N   is  the  number  of  time  points  (power  of  2), 

•  Y(k)     is  the  FFT  of  the  signal  y(t)  and 

•  S(k)     is  the  periodogram  spectral  estimate  of  y(t) . 

This  procedure  was  developed  through  the  computer  program 
SPECTRUM. M  in  MATLAB.  The  program  listing  appears  in 
Appendix  D. 


B.   SYNTHETIC  SIGNAL  TESTING  RESULTS 

This  procedure  shows  the  weaknesses  of  the  Cadzow-Solomon 
algorithm  by  using  no  bias  compensation.  Cases  were  examined 
for  the  three  different  SNR's  of  30,  20  and  10  dB. 
Overestimation  varied  from  2:1  up  to  5:1,  e.g.,  for  20  true 
poles  and  60  asking  poles,  the  ratio  is  3:1.  The  interval 
processed  was  composed  of  2  00  points,  starting  where  the 
excitation  began. 

Figures  11-23  illustrates  the  results  of  this  effort. 

Although  the  generated  SNR's  were  30,  2  0  and  10  dB,  the 

resulting  SNR's  over  the  processed  window  were  2  dB  higher. 

Observations  made  led  to  the  following  factors,  thus  improving 

results   as  previously  obtained   from  the  Cadzow-Solomon 

algorithm. 

•  The  most  significant  factor  is  to  select  the  excitation 
starting  point  as  the  point  to  start  the  algorithm 
processing. 
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•  The  second  significant  factor  is  to  select  the  system 
order  or  number  of  asking  poles.  In  the  case  of  20  true 
poles,  60  asking  poles  gives  the  most  accurate  results  for 
all  SNR's  employed. 

•  A  third  significant  factor  is  to  determine  the  number  of 
unknown  zeros.  For  synthetic  data,  the  position  of  the 
low  frequency  poles  was  obtained  with  better  accuracy  when 
this  number  was  equal  to  the  number  of  asking  poles  than 
with  the  case  where  the  number  of  asking  zeros  was  equal 
to  the  number  of  true  poles. 

Figure  23  illustrates  the  case  of  K+l=20  (as  compared  to 

Figure  19  for  K+l=60) ,  where  (K+l)  represents  the  number  of 

asking  zeros.    In  the  case  where  K+l=20,  more  poles  were 

obtained  within  the  unit  circle  (as  compared  to  the  case  where 

K+l=60)  at  the  frequencies  of  the  true  poles,  but  with  a 

different  damping  factor.  The  experimental  data  indicates  the 

following  result: 

•  There  is  a  bias  of  all  the  poles  towards  the  unit  circle 
that  results  in  the  loss  of  some  poles.  The  bias  is  so 
large  that  it  forces  the  poles  to  appear  on  the  other  side 
of  the  unit  circle  as  noise  poles. 

By  examining  the  spectra  of  the  synthetic  signal  at  different 

noise  levels,  as  illustrated  in  Figures  15-16,  it  may  be 

noticed  that  the  high  frequency  poles  cannot  be  completely 

separated  from  the  white  noise  spectrum  when  increasing  the 

noise  level .    This  is  due  to  the  fact  that  these  high 

frequency  poles  carry  very  little  energy. 

The  set  of  equations  in  matrix  form  developed  in  Equation 

(29)  have  been  developed  in  the  computer  program  CADS0L1.M  in 

MATLAB.   The  program  listing  appears  in  Appendix  E. 
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Figure   11         Double-Gaussian  Smoothing  Function 
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Figure   12         Spectrum  of  the  Double-Gaussian  Function 
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Figure  14    Synthetic  Signal  with  SNR=10  dB 
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Figure   15         Spectrum  of  the  synthetic  signal  without  noise 
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Figure   16         Spectrum  of  the  synthetic   signal  with  SNR=10   dB 
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Figure   19  Poles  in  s-plane,   SNR=22  dB,    for  synthetic  signal 
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Figure  20         Poles  in  z -plane,   SNR=22  dB,   for  synthetic  signal 
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Figure   21         Poles  in  s-plane,   SNR=12  dB,    for  synthetic  signal 
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C.   THIN  WIRE  SIGNAL  TESTING 

The  Cadzow-Solomon  algorithm  was  also  tested  using  thin 
wire  measured  scattering  data.  The  results  have  been  compared 
with  the  results  obtained  using  time  domain  integral  equation 
(TDIE)  thin  wire  data.  The  poles  obtained  by  processing  the 
TDIE  computed  data  were  assumed  to  be  correct,  even  though  the 
program  that  computes  the  currents  on  the  thin  wire  does  not 
take  into  account  the  capacitance  at  the  ends  of  the  wire. 
1.  Thin  Wire  Integral  Equation  Computed  Data 

For  the  thin  wire,  the  natural  resonances  may  be 
determined  by  solving  the  integral  equation  that  describesthe 
current  flowing  on  the  wire.  The  result  of  this  type  of 
simulation  is  the  response  of  the  wire  to  a  specified 
excitation  field.  Morgan  [Ref.  14]  developed  a  time- 
domain  thin  wire  integral  equation  computer  routine,  based  on 
the  formulations  of  Sayre  and  Harrington  [Ref.  15]. 
The  wire  used  for  this  simulation  had  a  length  of  0.1  meter 
and  a  radius  of  1.18  mm.  The  back-scattering  response  was 
computed  at  three  different  incident  aspects,  ranging  from  30 
to  90  degrees,  with  90  degrees  representing  a  broadside 
aspect.  The  excitation  waveform  used  was  the  same  double 
Gaussian  as  that  used  with  the  synthetic  data,  as  illustrated 
in  Figure  11.  Figure  30  illustrates  the  position  pole 
estimates  obtained  using  Cadzow-Solomon  algorithm  and  shows 
the  aspect  independence  of  the  poles  for  three  cases.   The 
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five  low  frequency  poles  appeared  exactly  at  the  same  position 
for  all  aspects.  It  should  be  noted  that  these  results,  as 
illustrated  in  Figure  30,  appeared  to  be  exactly  the  same, 
irrespective  of  any  variations  in  the  parameters  (No.  of 
points  -  No.  of  asking  poles  -  No.  of  asking  zeros)  used  in 
processing  the  signal.  Note  that  for  broadside  illumination, 
only  the  odd-numbered  poles  appear.  Because  of  the  physical 
symmetry  of  both  the  thin  wire  and  incident  field,  the  even- 
numbered  modes  are  prevented  from  being  excited  by  the 
incident  field.  With  illumination  at  45  degrees  aspect,  a 
spectrum  with  adequate  energy  only  within  the  bandwidth  from 
1  to  8  GHz  was  produced,  as  expected.  At  higher  frequencies, 
backscattering  is  suppressed  because  most  of  the  energy  is 
reradiated  near  to  the  specular  scattering  directions. 
Therefore,  only  the  five  low  frequency  poles  are  accurately 
obtained  for  this  case. 

The  TDIE  program  generated  a  response  at  30  degrees 
aspect  consisting  of  255  points  over  5  nsec,  resulting  in  a 
sampling  frequency  of  50.8  GHz.  In  the  case  of  45  and  90 
degrees  aspect  consisting  of  240  points  over  5  nsec,  a 
sampling  frequency  of  47.8  GHz  resulted.  These  sampling 
frequencies  differ  from  the  sampling  frequency  of  51  GHz  used 
in  the  measured  data. 

When  processing  the  TDIE  thin  wire  data,  points  from 
160  to  200  were  used  to  run  the  algorithm  starting  at  the 
point  where  the  excitation  starts.  The  number  of  asking  poles 

53 


ranged  from  28  (for  the  90  degrees  aspect)  to  40  (for  the  30 
degrees  aspect) .  The  number  of  asking  zeros  or  feed-forward 
order  of  the  system  was  either  the  same  as  the  number  of 
asking  poles  or,  was  calculated  by  determining  the  length  of 
early-time.   This  early-time  was  computed  from  the  formula 

t  =  —  -(1+cosa)  (49) 

e  c 

where , 

•  L   is  the  length  of  the  wire  and 

•  a  is  the  aspect  angle  from  end-on  orientation. 

This  value  of  time  was  then  converted  to  the  appropriate 
number  of  time  points,  as  based  on  the  sampling  time  interval 
of 


t. 


+1 


(50) 


nb=  integer 

where, 

•  T  is  the  sampling  interval,  and 

•  nb   is  the  feed-forward  order  of  the  system. 

This  value  of  nb   is  the  minimum  number  of  asking  poles  that 

may  be  used,  presenting  the  number  of  delays  in  the  z- 
transform  for  the  early-time  of  the  system. 

Another  consideration  in  the  processing  of  these  data 
was  the  scaling  of  the  excitation  waveform  and  its  position 
with  respect  to  the  computed  response.  Although  the  data  were 
generated  using  a  double  Gaussian  waveform  with  a  1  Volt  peak 
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amplitude,  the  excitation  waveform  had  approximately  the  same 
peak  amplitude  as  the  response  waveform.  Such  scaling  does 
not  change  the  frequency  contents  of  the  exciting  waveform  and 
gives  better  results  in  the  pole  extraction  due  to 
minimization  of  some  of  the  effects  of  ill-conditioning  in  the 
data  matrix. 

To  position  the  driving  waveform  with  respect  to  the 
response  waveform,  the  time  difference  between  the  excitation 
of  the  first  and  last  point  of  the  wire  may  be  computed.  This 
time  interval  is  represented  as 

tdelay=—' COSa,  (51) 

where  a  is  the  aspect  angle.  The  maximum  absolute  value  of 
the  response  occurs  when  the  last  point  of  the  wire  is  being 
excited.  The  information  derived  from  the  time  interval, 
tdelay,    and  the  maximum  absolute  value  of  the  response  allows 

the  excitation  waveform  to  be  positioned  with  the  first  point 
of  the  wire. 

2 .  Thin  Wire  Measured  Data 

Measurements  were  performed  in  the  Transient 
Electromagnetic  Scattering  Lab  (TESL)  to  test  the  algorithm  by 
using  a  thin  wire  with  the  same  dimensions  as  the  one  used  for 
the  previous  simulation.  A  detailed  explanation  of  the 
procedure  and  the  techniques  used  for  measuring  the  scattering 
response  from  the  thin  wire,  as  well  as  other  scale  model 
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targets  may  be  found  in  Morgan  and  McDaniel   [Ref.  16] 
and  Bresani  [Ref.  17]. 

One  measurement  was  available  for  each  of  the  aspects 
at  30,  45  and  90  degrees.  The  scattering  response  waveforms 
obtained  from  the  measurements  are  illustrated  in  Figures  24, 
26  and  28.  These  waveforms  were  compared  with  the  calculated 
waveforms  obtained  from  the  TDIE  program.  It  can  be  seen  from 
the  figures  that  the  natural  modes  between  the  calculated  and 
the  measured  waveforms  do  not  have  exactly  the  same 
frequencies. 

The  spectrum  of  the  measured  response  waveforms  were 
also  obtained,  giving  a  distribution  of  energy  within  the 
bandwidth  1-12  GHz  on  eight  frequencies,  as  shown  in  Table 
III: 


Table  III 

THIN  WIRE 


DISTRIBUTION  OF  ENERGY  WITHIN  THE  SPECTRUM  OF 


Aspect 
[degrees] 

Frequencies 
[GHz]  with 
most  energy 

Frequencies 
[GHz]  with 
less  energy 

Frequencies 
[GHz]  with 
no  energy 

30 

3,4.4,5.8, 
7.2 

1.6,8.6,10.2 

11.6 

45 

3,4.4 

1.6,5.6 

7.2,8.6,10.2, 
11.6 

90 

1.6,4.4 

7.2 

3,5.6,8.6, 
10.2,11.6 
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Figures  25 ,  27  and  29  illustrate  the  spectrum  of  the  measured 
thin  wire  scattering  response  for  each  aspect. 

When  the  measured  thin  wire  response  was  processed, 
the  points  used  to  run  the  algorithm  were  from  180  to  200, 
starting  at  the  initial  excitation  point,  while  the  number  of 
asking  poles  ranged  from  30  to  60.  The  best  results  were 
obtained  when  the  number  of  asking  poles  approached  40.  The 
feed-forward  order  of  the  system  was  either  the  same  as  the 
number  of  asking  poles,  or  was  calculated  by  determining  the 
early-time  length.  The  scaling  of  the  driving  waveform  and 
the  positioning  of  the  excitation  with  respect  to  the  response 
was  computed,  as  previously  described.  Figure  3  0  illustrates 
the  extraction  of  the  poles  from  the  TDIE  response,  while 
Figure  31  illustrates  the  extraction  of  the  poles  from  the 
measured  waveform  for  the  combined  aspects.  Figures  32-37 
illustrate  the  extraction  of  the  poles  from  the  measured 
waveform  for  each  aspect  separately.  The  pole  results 
obtained  using  both  the  bias  compensation  method  developed  in 
this  thesis  and  Cadzow-Solomon's  bias  compensation  method  have 
been  plotted  in  Figures  32-37  along  with  the  poles  obtained 
without  any  bias  compensation  method,  so  that  the  results  from 
the  three  types  of  methods  may  be  compared. 
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Figure   24  Thin  wire  TDIE  &  Measured  Response,    3  0  deg  aspect 
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Figure  25    Spectrum  of  thin  wire  measured  response,  3  0  deg 
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Figure   2  6         Thin  wire  TDIE  &  Measured  Response,   45  deg  aspect 
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solid  :  computed  by  TDIE  dashed  :  measured        dotted  :  excitation 

Figure   28  Thin  wire  TDIE  &  Measured  Response,    90  deg  aspect 
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Figure  35    Poles  at  z-plane,  30  deg  aspect 
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IV.  POLES  FROM  SCALE  MODEL  AIRCRAFT 

Scattering  data  for  four  different  scale  model  aircraft 
targets  were  processed  without  bias  compensation,  using  the 
Cadzow-Solomon  algorithm.  Poles  were  extracted  from  the 
measured  scattering  responses  of  the  aircraft  targets  to 
double-Gaussian  electromagnetic  excitation  incident  at  0,  30 , 
90,  and  180  degrees.  The  0  degree  aspect  represents  a  nose-on 
measurement,  while  the  90  degree  represents  a  broadside 
measurement.  Signatures  are  shown  in  Figures  38-45  for  targets 
1  and  2 . 

The  bandwidth  of  the  TESL,  1-12  GHz,  was  matched  to  the 
scaling  factor  of  the  model  aircraft  targets.  This  scaling 
factor  was  1/72  for  all  of  the  model  aircraft  targets  used. 
Table  IV  contains  the  full  scale,  significant  dimensions  of 
these  aircraft  targets.  The  aircraft  data  was  collected  by 
Bresani  [Ref.  17]  at  the  four  incident  angles,  as  previously 
described. 


Table  IV     SCALED  DIMENSIONS 

OF 

AIRCRAFT 

TARGETS 

Target  number            1 

2 

3        4 

Overall  length  (cm)      26.5 

22.5 

23.6     26.2 

Wingspan  (cm)            19.8 

19.0 

16.4     16.0 
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A.  DATA  PROCESSING 

The  model  aircraft  scattering  data  was  processed  using  a 
number  of  time  points  ranging  from  200  to  240.  The  algorithm 
processing  was  started  at  the  point  of  initial  excitation. 

The  positioning  of  the  excitation  waveform  was  set  up 
manually  with  respect  to  the  response  waveform  from  observa- 
tions of  the  scattering  data.  Manual  positioning  for  the 
driving  waveform  was  required  as  no  obvious  condition  existed 
for  the  model  aircraft,  as  was  the  case  for  the  thin-wire. 

The  number  of  asking  poles  was  set  to  60,  as  the  value  of 
60  was  found  from  previous  experimentation  to  represent  the 
most  suitable  results.  The  number  of  asking  zeros  was  the 
same  as  the  number  of  asking  poles.  This  value  was  based  on 
the  fact  that  the  computed  early-time  for  the  aircraft  target 
is  usually  about  2L/c.  Conversion  of  this  early-time  to  the 
appropriate  number  of  time  points  (equal  to  the  number  of 
asking  zeros)  provides  a  number  larger  than  the  number  of 
asking  poles.  However,  the  number  of  asking  zeros  may  not  be 
larger  than  the  number  of  asking  poles,  the  two  values  are  set 
to  be  equal. 

B.  RESULTS  FROM  EXTRACTING  THE  POLES 

The  model  aircraft  scattering  data  showed  that  the  poles 
were  less  likely  to  group  than  the  thin-wire  data.  The  highly 
complex  nature  of  the  aircraft  scatterer  combined  with  the 
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small  number  of  significant  data  points  (about  150)  became  a 
difficult  problem  for  the  algorithm  to  solve. 

Figures  46-47  illustrate  the  spectrum  of  the  data  record 
for  aircraft  targets  1  and  2,  for  different  incident  angles. 

Figures  48-51  illustrate  the  extraction  of  the  poles  from 
the  same  aircraft  target  for  the  combined  aspects.  To  obtain 
an  initial  indication  of  the  possibility  of  target  identifi- 
cation by  pole  extraction,  nose-on  measurement  results  have 
been  compared  in  Figures  52-53. 

The  thin  wire  data  (Chapter  III)  indicated  that  the 
principle  poles  provided  well  defined  clusters  for  all  the 
incident  angles  examined,  despite  the  use  of  a  wide  range  of 
processing  parameters.  However,  the  scale  model  aircraft 
targets  did  not  provide  any  well  defined  clusters  under  the 
limited  attempts  at  pole  estimation  conducted  here.  The  prin- 
ciple reason  for  this  difference  appears  to  be  that  the  scale 
models  have  more  complicated  pole  patterns  than  the  thin  wire 
targets  and  insufficient  time  was  available  to  explore  ideas 
for  optimizing  the  performance  of  the  estimation  algorithm. 

Processing  the  experimental  scattering  data  for  aircraft 
models,  as  performed  herein,  is  only  an  initial  attempt  to 
demonstrate  resonance  estimation  for  real-world  targets 
embodying  complex  configurations.  Processing  for  this 
aircraft  data  was  conducted  for  only  one  week  at  the  end  of 
the  thesis  research.  It  is  recommended  that  this  same  data  be 
used  for  a  follow-on  thesis. 
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Figure   38 
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Figure  39    Measured  response,  model  aircraft  1,   30  deg 
aspect 
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solidrresponse  dashrexcitation 

Figure  40    Measured  response,  model  aircraft  1,   90  deg 
aspect 
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Figure  41    Measured  response,  model  aircraft  1,  180  deg 
aspect 
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Figure   42         Measured  response,  model  aircraft  2,    0  deg  aspect 
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Figure  44    Measured  response,  model  aircraft  2,   90  deg 
aspect 
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Figure  45    Measured  response,  model  aircraft  2,  180  deg 
aspect 
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Figure  46    Spectrum  of  model  aircraft  1  measured  response, 
all  aspects 
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Figure  47    Spectrum  of  model  aircraft  2  measured  response, 
all  aspects 
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Figure  48    Poles  extracted  from  model  aircraft  1  measured 
response  in  s-plane 


77 


—                                          X 

+? 

+ 

- 

-  •  .. 

o 

+ 

0 

M 

* 

+  :0  deg 

• 

<9. 

- 

o:30  deg 

M 

- 

x90  deg 

X 

M80  deg 

• 

- 

M 

* 

0.9- 


0.8 


0.7- 


0.6 


0.5 


0.4- 


0.3- 


0.2- 


0.1- 


Figure  49    Poles  extracted  from  model  aircraft  1  measured 
response  in  z-plane 
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Figure  50    Poles  extracted  from  model  aircraft  2  measured 
response  in  s-plane 
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Figure  51    Poles  extracted  from  model  aircraft  2  measured 
response  in  z-plane 
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Figure  52    Poles  extracted  from  four  model  aircraft  in  s- 
plane,  nose-on 
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Figure  53    Poles  extracted  from  four  model  aircraft  in  z- 
plane,  nose-on 
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V.  SUMMARY  AND  CONCLUSIONS 


A.   SUMMARY 

This  thesis  has  attempted  to  demonstrate  radar  target 
identification  by  building  on  the  earlier  work  of  Norton  [Ref . 
11],  Larison  [Ref.  12]  and  Murphy  [Ref.  13].  The  largest 
portion  of  the  work  consisted  of  testing  the  Cadzow-Solomon 
algorithm  using  synthetic  and  thin  wire  measurement  data. 

Chapter  I  introduced  the  use  of  the  Kumaresan-Tufts  and 
the  Cadzow-Solomon  algorithms  to  locate  target  poles  for  a 
Non-Cooperative  Target  Recognition  system.  The  resonance- 
based  radar  target  identification  problem  was  discussed  and 
the  transient  electromagnetic  scattering  described. 

Chapter  II  consisted  of  two  parts.  The  first  part 
discussed  early  methods  used  to  solve  target  identification 
problems.  This  discussion  included  a  description  of  Prony's 
method  and  also,  Singular  Value  Decomposition  on  which  the 
Kumaresan-Tufts  and  Cadzow-Solomon  algorithms  were  developed. 
The  Kumaresan-Tufts  algorithm  was  developed  in  detail, 
including  Norton's  bias  compensation  method.  The  performance 
of  the  Kumaresan-Tufts  algorithm  was  also  demonstrated,  as 
illustrated  in  Figures  3  through  10.  The  second  part 
described  the  Cadzow-Solomon  algorithm.  Two  bias  compensation 
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methods  were  included:  that  of  Cadzow-Solomon  and  that  of  the 
author. 

Chapter  III  demonstrated  the  Cadzow-Solomon  algorithm 
testing  in  two  phases.  The  first  phase  of  testing  was 
performed  with  synthetic  data,  while  the  second  phase  was 
performed  with  thin  wire  measurement  data.  The  thin  wire  data 
testing  phase  attempted  to  evaluate  conditions  appearing  from 
a  real  target  response.  Results  of  the  synthetic  signal  model 
testing  are  illustrated  in  Figures  13  through  23.  The  results 
of  the  thin  wire  scattering  signature  testing  are  illustrated 
in  Figures  24  through  37. 

Chapter  IV  considered  an  initial  attempt  at  extraction  of 
poles  from  scale  model  aircraft  measurements  obtained  in  the 
NPS  Transient  Electromagnetic  Scattering  Lab.  Measurement 
data  was  processed  using  the  Cadzow-Solomon  algorithm  and  the 
unsuccessful  results  are  illustrated  in  Figures  38  through  53. 

Testing  of  both  the  Kumaresan-Tufts  and  the  Cadzow-Solomon 
algorithms  was  performed  using  MATLAB  programs.  A  sequence  of 
programs  was  written  to  complete  the  demonstration  of  the 
target  identification  problem.  Appendices  A  to  E  present  only 
a  part  of  this  sequence  of  programs,  including  the  theoretical 
models  of  the  Cadzow-Solomon  algorithm. 
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B.   CONCLUSIONS 

Both  the  Kumaresan-Tufts  and  the  Cadzow-Solomon  algorithms 
can  effectively  extract  poles  from  the  scattering  response  of 
simple  high-Q  targets  such  as  thin  wires.  As  the  late-time 
portion  of  the  response  signal  is  weak  (with  low  SNR)  and  the 
Cadzow-Solomon  algorithm  has  the  ability  to  use  the  early-time 
portion,  this  thesis  concentrated  mainly  on  the  use  of  the 
Cadzow-Solomon  algorithm. 

The  Cadzow-Solomon  algorithm  extracts  poles  of 
synthetically  generated  data,  integral  equation  computed  data, 
and  thin  wire  measured  scattering  data  with  excellent  results. 
A  signal-to-noise  ratio  above  10  dB  is  required,  depending  on 
the  damping  rate  of  the  data.  The  system  order  is 
intentionally  overestimated.  The  excess  order  provides  the 
flexibility  to  model  the  noise  and  improves  the  estimation  of 
parameters  of  exponentially  damped  sinusoidal  signals  in 
noise.  The  Singular  Value  Decomposition  method  alleviates 
severe  ill-conditioning  of  the  data  matrix.  Backward 
prediction  and  SVD  are  used  to  separate  the  computed  signal 
poles  from  spuriously  computed  noise  poles  introduced  by  the 
overestimated  system  order. 

The  most  critical  parameters  required  for  the  successful 
thin-wire  processing  were  the  selection  of  the  appropriate 
starting  point  to  begin  the  data  processing  and  the 
appropriate  system  order.  The  best  results  were  obtained  with 
the  starting  point  of  the  data  processing  set  at  the 
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excitation  starting  point,  and  with  the  processing  system 
order  at  about  3  times  the  true  system  order. 

The  existence  of  noise  in  the  target's  response  produces 
a  bias  in  the  positioning  of  the  extracted  poles.  Thus, 
several  bias  compensation  schemes  may  be  developed. 

Given  a  very  limited  initial  effort,  the  Cadzow-Solomon 
algorithm  was  unable  to  extract  poles  of  scale  model  aircraft 
scattering  data  with  adequate  accuracy.  It  is  conjectured 
that  the  data  points  in  the  natural  mode  information  response 
are  too  few  for  the  algorithm  to  model  both  the  target  poles 
and  the  noise  poles. 
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APPENDIX  A.  ARMA  COEFFICIENT  AND  RECURSIVE  SIGNAL  GENERATOR 


PROGRAM  LISTING 

%  ARMAGEN1.M 

% 

%   Program  generates  a(k),  b(k)  coefficients  of 

% 

%  b(0)+b(l)*zA-l+. . .+b(q)*zA-q 

%       H(z)  =  

%  1  +a(l)*zA-l+. . .+a(p)*zA-p 

% 

%    given  the  s-plane  poles,  residues,  number  of  time  points 

%   and  time  window. 

%    Programed  by  Gregory  Lazarakos,  5  Mar  1991. 

% 

!del  temp. mat 

case=  * synt ' ; 

dt=tO./(notp-l) ; 

i=sqrt(-l) ; 

dm=exp( (sigma+i* omega) *dt) ; 

alpha=ampl. *exp (i*phase) ; 

dm=dm' ; 

alpha=alpha' ; 

[ 'Please  wait. . . ' ] 

[b,a]=residue(alpha,dm, 0) ; 

A=real (a) ; 

B=real(b) ; 

%    Program  generates  time  response  of  an  ARMA  system 

%   via  the  equation 

%  N  L 

%   y(n)  =  SUM  A(k) *y (n-k+1)  +  SUM  B(k) *x(n-k+l)  ,forn=l:notp 

%  k=2  k=l 

%   given  A(k) ,  B(k)  and  input  excitation  record  x(n) . 

%    Programed  by  Gregory  Lazarakos,  10  Mar  1991. 

N=size(A) ; 

N=N(2) ; 

L=size(B) ; 

L=L(2); 

for  j=2:N 

A(j)=-A(j); 

end 
% 

[•Excitation  is  double  Gaussian'] 

excgen 
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% 


[ • Please  wait. . . ■ ] 


ye=zeros(l/notp) ; 
ye(l)=B(l)*x(l); 
energy=(ye(l) ) . A2 ; 
for  n=2 : notp 

ye (n) =0.0; 

Ln= [ n ; L] ; 

Lmax=min(Ln) ; 

for  k=l:Lmax 

ye(n)=ye(n)+B(k)*x(n-k+l) ; 

end 

Nn=[n;N] ; 

Nmax=min(Nn) ; 

for  k=2:Nmax 

ye(n)=ye(n)+A(k)*ye(n-k+l) ; 

end 

energy=energy+ (ye (n) ) . A2 ; 
end 

rms=sqrt (energy) ; 
ic=zeros (l,notp) ; 
for  n=l:notp 

ic(n)=ye(n) ./rms; 
end 

axis([0  notp  -0.4  0.8]) 
plot (ic) 
title ([ 'Generated   signal   w/o   noise,   medium   Q,   by 
' fint2str(nop) , •  poles']); 
xlabel('time  points'); 
ylabel ( ' amplitude  rms ' ) ; 
grid  ; 
pause 

hard=input  ( '  Do  you  want  a  hardcopy  for  the  plot?  [n]  : 
\'s'); 

if  hard=='y' 

plotfile=input( 'Enter  filename  :  ','s'); 
eval ( [ ' meta  ' , plotf ile ] ) 
end 


% 
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APPENDIX  B.  DOUBLE  GAUSSIAN  SMOOTHING  FUNCTION  GENERATOR 


PROGRAM  LISTING 

%  EXCGEN.M 

% 

%   Program  generates  excitation  record  x(n)  for  a  response 
%    from  TESL  data  and  synthetic  data, 
%   via  the  equation  of  a  double  Gaussian  waveform 
% 

%   x(n)  =  Al*exp(-al*tA2)  -  A2*exp(-a2*tA2) 
% 

%   given  the  time  window. 

%    Program  uses  input  values  of  the  10%  height  of  the  low 
%    and  high  frequency  ends  of  the  double  Gaussian  frequency 
%   response  to  determine  the  pulse  widths  in  the  time  domain. 
%    Programed  by  Gregory  Lazarakos,  8  Apr  1991. 
%    Last  Revision  August  6,  1991. 
if  case=='meas' 

if  filename (1)=='F' 
LENGTH=0.25; 

elseif  filename(siz-l:siz)=='sp' 
LENGTH=eval ( filename (3:4)); 
LENGTH=LENGTH/100 ; 

else 

LENGTH=0 . 1 ; 

end 
end 

npts=notp ; 

tmin=0. 0; 

tmax=t0 ; 

DT=dt ; 

Tl=0.147; 

T2=0.314; 

thr=input( 'Enter  threshold  value   for  time   in   nsec 
[1.25]: •) ; 

if  isempty(thr) 
thr=1.25; 

end 

al=(4.0.*log(10)) ./(T1.A2) ; 

a2=(4.0.*log(10)) ./(T2.A2)  ; 

Al=sqrt(al) ./ (sqrt (al) -sqrt (a2) ) ; 

A2=A1-1.0; 
if  case=='synt' 

point=input( 'Which  to  be  the  point,   the  excitation 
starts? • ) ; 

x=zeros(l,npts) ; 

for  n=l:npts 
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t=(n-ll-point) .*DT; 
if  t<=thr 

el=Al.*exp(-al.*t.A2)  ; 
e2=A2.*exp(-a2.*t.A2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 
end 
end 
elseif  case==lmeas' 

asprad=aspect* (pi/180)  ; 
top=input ( *  Do  you  want  to  position  the  excitation,  auto  or 
manual  ?  {a,m}  :  *  , 's'); 
if  top=='a' 

[maxic , indexl ] =max ( ic) ; 
[minic, index2]=min(ic) ; 
x=zeros(l,npts) ; 

tdel=cos(asprad)*(2*LENGTH/3e8)  ; 
ndel=round(tdel/(dt*le-9) ) ; 
if  abs (maxic) >=abs (minic) 
point=indexl-ll-ndel ; 
for  n=l:npts 

t=(n-indexl+ndel) .*DT; 
if  abs(t)<=10*DT 

el=Al.*exp(-al.*t.A2) ; 
e2=A2.*exp(-a2.*t.A2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 
end 
end 

x=x*0.6; 
else 

point=index2-ll-ndel ; 
for  n=l:npts 

t=(n-index2+ndel) .*DT; 
if  abs(t)<=10*DT 

el=Al.*exp(-al.*t.A2) ; 
e2=A2.*exp(-a2.*t.A2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 
end 
end 

x=x*0.6; 
end 
end 
if  top=='m' 

point=input( 'Which  to  be  the  point,  the  excitation 
starts  ?  • ) ; 

x=zeros(l,npts) ; 
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for  n=l:npts 

t=(n-ll-point) .*DT; 
if  abs(t)<=10*DT 

el=Al.*exp(-al.*t.A2) ; 
e2=A2.*exp(-a2.*t.A2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 
end 
end 

x=x* . 6 ; 
end 

ltst= ( 1+cos (asprad) ) * (LENGTH/3e8 ) ; 
nltst=round(ltst/(dt*le-9) )  ; 
ltimp=point+22+nltst; 
else 

exit 
end 

if  case== • synt ' 

axis([0  notp  -0.4  1]) 
plot(x) 
title ( 'Double-Gaussian  Smoothing  Function'); 
xlabel ( • time  points ' ) ; 
ylabel ( 'amplitude  rms'); 
elseif  case=='meas' 

plot (l:npts,x, ' — ' , l:notp, ic, ltimp, ic(ltimp) , ' * ' ) 
if  filename (1)=='F' 
title  ( [  •  Data  response  signal  from  aircraft '  ,  filename  (1:3)  ,  • , 
aspect=' ,num2str (aspect) , '  deg  ' ]) ; 

elseif  f ilename(siz-lrsiz) =='sp' 
title  ( [  '  Data  response  signal  from  sphere  '  ,  filename  (3:4)  ,  'cm, 
aspect=' ,num2str( aspect) , '  deg  ']); 
else 
title (['Data    response    signal    from    thin    wire, 
aspect=' ,num2str (aspect) , '  deg  ' ]) ; 
end 

xlabel ( 'solid=response        dash=excitation' ) ; 
ylabel ( 'amplitude  ' )  ; 

text  ( ltimp  ,0.1,  • late  time > '  >  ; 

if  filename (1: 2) =='et' 

text (100, 0.4, 'TDIE  response'); 
end 
grid  ; 
else 

exit 
end 

pause 
hard=input ( ' Do  you  want  a  hardcopy  for  the  plot? 
[n]:\'s'); 

if  hard=='y' 
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plotf ile=input( 'Enter  filename  :  * ,'s'); 
eval(['meta  •  , plotf ile]) 
end 
if  case==,meas' 

save  temp  ic  energy  x  point  dt  to  notp  filename  aspect  siz 
LENGTH 
end 
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APPENDIX  C.  SYNTHETIC  NOISE  GENERATOR 

PROGRAM  LISTING 

%  NOISEGEN.M 

% 

%        Synthetic  Noise  Generation 

%        Programmed  by  Greg  Lazarakos  20/2/91 

% 

SN=input('How  many  db  ?    (7,10,20,30)     [20]:'); 

if  isempty(SN) 
SN=20; 

end 

sn=SN./10; 

sn=10. Asn; 

en=0; 

for  k=l:notp 

en=en+ ( ic ( k) ) . A2 ; 

end 

ava=en . / (notp . *sn) ; 

ava=sqrt (ava) ; 

rand( 'normal ' ) 

w=rand ( ic) *ava ; 

ic=ic+w; 
% 

save  temp  ic  x  point  filename  w  ye  energy  dt  dm  notp  omega 
sigma  to 
% 

axis([0  notp  -0.4  0.8]) 

plot(ic) 
title  ( [  'Generated  signal  with  noise  S/N='  ,  int2str(SN)  ,  '  db ' ] ) ; 

xlabel('time  points'); 

y label ( ' amplitude  rms ' ) ; 

grid  ; 

pause 

hard=input  ( '  Do  you  want  a  hardcopy   for   the   plot? 
[n]:\'s'); 

if  hard=='y« 

plotfile=input( 'Enter  filename  :  ','s'); 
eval(['meta  ',plotfile]) 
end 
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APPENDIX  D.  SPECTRUM  ESTIMATION 

PROGRAM  LISTING 

%  SPECTRUM. M 

% 

%    Programed  by  Gregory  Lazarakos,  16  Apr  1991. 

%    Last  revision,  6  Aug  1991. 

df=l./((notp-l) .*dt) ; 

IC=fft(ic,notp) ; 

f=(0:notp-l)*df ; 

fmax=l./dt; 

axis([0  fmax  -3  3]) 

plot(f ,IC) 

title ( 'Fourier  transform  of  the  response  signal'); 

xlabel (' Frequency  in  GHz'); 

ylabel ( 'Amplitude ' ) ; 

grid; 

pause 

SIC=(abs(IC)) .  A2; 
SIC=SIC/notp; 

axis([0  20  0  max(SIC)]) 

plot(f,SIC) 

if  isempty(SN) 

title ( 'Spectrum  of  the  response  signal  •); 

else 

title ([ 'Spectrum    of    the    response    signal    with 
S/N=' ,int2str(SN) , •  db']); 

end 

xlabel (' Frequency  in  GHz'); 

ylabel ( ' Ampl itude ' ) ; 

grid; 

pause 

hard=input( 'Do  you  want  a  hardcopy  for  the  plot? 
[n]  :  \'s'); 

if  hard=='y' 

plotfile=input( 'Enter  filename  :  ','s'); 
eval(['meta  ',plotfile]) 
end 
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APPENDIX  E.  THE  CAD ZOW- SOLOMON  POLE  EXTRACTION  ALGORITHM 

PROGRAM  LISTING 

%  CADSOL1.M 

% 

%  Cadzow-Solomon's  Algorithm  for  Extracting  the  Poles 

%  and  Residues 

%  Version  1.0 

%  Programmed  by  Greg  Lazarakos  3/16/91 

% 

%  Set  the  first  sample  (kapa) 

if  point>0 

kapa=point; 

else 

kapa=l; 

end 
%  Set  the  number  of  samples  (CN) 

CN=input ( 'How  many  samples  to  process  ?  (>120)  [200]:'); 

if  isempty(CN) 
CN=2  00; 

end 
%  Compute  the  SNR  for  the  processed  time  window 

en=0; 

nen=0; 

for  k=kapa:CN 

en=en+(ic(k) ) . A2; 
nen=nen+(w(k) ) . A2 ; 

end 

SNR=en./nen; 

SNR=10.*logl0(SNR) ; 
%  Set  the  number  of  poles  (CL)  >  number  of  real  poles  (CM) 
%  CM  <=  CL  and  2*CL  <=  CN-CL 

CL=input( 'How  many  unknown  poles  ?  (>20)  [60]:  •); 

if  isempty(CL) 
CL=60; 

end 
%  Set  the  number  of  expected  zeros  (CK) 

CK=input ( 'How  many  expected  zeros  ?  [No.  poles]:  '); 

if  isempty(CK) 
CK=CL  " 

end 

[ 'Please  wait. . . ' ] 
%  Forming  the  matrix  YI[ (CN-CL) *CL]  from  the  data  ic(n) 
YI=zeros (CN-CL, CL) ; 

for  j=kapa+l:CN-CL+kapa 

YI(j-kapa, : )=ic(j : j+CL-1) ; 

end 
%  Forming  the  matrix  h[ (CN-CL)*!]  from  the  data  ic(n) 
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h=[  ]; 

h=ic(kapa:CN-CL+kapa-l) ; 
h=h  •  ; 
%  Forming  the  matrix  XI [ (CN-CL) *CK]  from  the  data  x(n) 
XI=zeros (CN-CL, CK) ; 
for  j=l: CN-CL 

XI(j,:)=x(j:j+CK-l); 
end 
%  Unify  matrices  YI  and  XI  as  I=[YI|XI] 

I=[YI  XI] ; 
%  Set  the  tolerance 

tol=0. 000001; 
%  Find  the  vector  of  backward  prediction  coefficients 

beta=pinv(I,tol) *[h] ; 
%   Set  the  coefficients  of  the  prediction-error  filter 
polynomial 

ca=zeros(l/CL+l) ; 
ca(l)=l; 
for  j=l:CL 

ca(j+l)=-beta(j)  ; 
end 
%  Set  the  coefficients  of  the  polynomial  B(z) 
cb=zeros(l,CK) ; 
for  j=l:CK 

cb(j)=beta(j+CL) ; 
end 
%  Find  the  residues  and  the  poles 

[resid/d/gk]=residue(cb/ca) ; 
%  Find  the  signal  poles 
s=log(d)/dt; 
k=0; 
1=0; 
for  j=l:CL 

if  real(s(j) )>0, 
k=k+l; 

polesl(k)=-s( j)  ; 
residl (k) =resid ( j ) ; 
else 

1=1+1 ; 

nopoll(l)=-s(j) ; 
noresl(l)=resid( j) ; 
end 
end 

dml=exp(polesl*dt) ; 
dm2=exp(nopoll*dt)  ; 
alphal=exp(residl*dt) ; 
CM=k; 

[  'The  signal  has  the  following  '  /num2str(CM)  ,  '  real  poles.  '  ] 
polesl 
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